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ABSTRACT 

The nature of angular momentum transport in the boundary layers of accretion disks 
has been one of the central and long-standing issues of accretion disk theory. In this 
work we demonstrate that acoustic waves excited by supersonic shear in the boundary 
layer serve as an efficient mechanism of mass, momentum and energy transport at the 
interface between the disk and the accreting object. We develop the theory of angular 
momentum transport by acoustic modes in the boundary layer, and support our findings 
with 3D hydrodynamical simulations, using an isothermal equation of state. Our first 
major result is the identification of three types of global modes in the boundary layer. 
We derive dispersion relations for each of these modes that accurately capture the 
pattern speeds observed in simulations to within a few percent. Second, we show that 
angular momentum transport in the boundary layer is intrinsically nonlocal, and is 
driven by radiation of angular momentum away from the boundary layer into both the 
star and the disk. The picture of angular momentum transport in the boundary layer 
by waves that can travel large distances before dissipating and redistributing angular 
momentum and energy to the disk and star is incompatible with the conventional notion 
of local transport by turbulent stresses. Our results have important implications for 
semianalytical models that describe the spectral emission from boundary layers. 

Subject headings: accretion, accretion disks - hydrodynamics — waves - instabilities 

1. Introduction 

Understanding the structure and dynamics of the boundary layer (BL) is one of the outstanding 
theoretical problems in accretion disk theory. We define the BL as the thin transitional region where 
the disk attaches to the star, gas slows down from the orbital to the stellar rotation velocity, and 
dVtjdw > is the angular velocity of the fluid and w is the cylindrical radius). In this work, we 
do not differentiate between the different types of accreting objects (a neutron star, a white dwarf, 
or a young star), as long as they have a surface. This assumes that the compact object has a weak 
enough magnetic field that the disk extends all the way down to the surface of the star without 
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being disrupted and channeled to the magnetic poles. The criteria for this to occur is that the 



magnetic pressure should be less than the ram pressure of the gas at the surface of the star (Ghosh 
fc Lamb|[l978| ). 



Within the BL, material slows down from the Keplerian rotational velocity to the rotational 
velocity inside the star. Since the radial width of the BL is significantly less than that of the star, 
this implies that the BL is a region of intense shear. In fact, the rate of energy dissipation inside 
the BL is expected to be as high as for the rest of the disk combined ( Kluzniak||1987 ), resulting in 
intense heating of the gas. 

An open question in BL theory is the mechanism of angular momentum transport. Due to a 
rising rotation profile, the BL is linearly stable to the magnetorotational instability (MRI) which 
is responsible for accretion in the disk (Balbus &; Hawley 1991). Indeed, Pessah fc Chan| Q2012 ) 
has shown that if dVtjdw > 0, then the stresses due to a sheared magnetic field are inefficient 
at transporting angular momentum and oscillate about zero. Thus, some other mechanism must 
be responsible for angular momentum transport within the BL. As we have previously shown in 



Belyaev et al. (2012), a likely candidate for this mechanism is excitation of acoustic waves in the BL 



accompanied by radiation of angular momentum away from the layer. This acoustic radiation arises 
due to sonic instabilities in the BL, which are a type of shear instability akin to the Papaloizou- 



flows ( |Glatzel|[T988| |Belyaev fc Rankov||2012l ) . 



Pringle instability (Papaloizou & Pringle 1984; Narayan et al. 1987) that operate in supersonic 



The focus of our present work will be the detailed description of acoustic modes in the BL 
and their relation to angular momentum transport and mass. We will study these modes both 
theoretically by solving for their dispersion relation and computationally with 2D, 3D unstratifled, 
and 3D stratified hydro simulations of the BL with an isothermal equation of state. 

Understanding the physics of acoustic modes is an important step towards understanding the 
temporal variability in the BL and constructing semianalytic models of BL spectra. Currently, 
virtually all semianalytic models of the BL assume that angular momentum transport is mediated 
by local turbulent stresses and adopt an effective ^-prescription based on Navier-Stokes viscosity 
( Popham fc Narayan||1995 ; Inogamov & Sunyaev|1999l|Piro fc Bildsten|2004aj ) . We argue, however, 
that angular momentum in the BL is transported by global waves, rather than local turbulence. 
These waves can potentially travel long distances before dissipating, meaning that angular momen- 
tum and energy transport is inherently nonlocal and cannot be well-described by any model based 
on a local anomalous viscosity. 

The paper is organized as follows. In ^2j we describe the model setup and numerical methods 
that we use for studying the BL, and in ^3] we provide a general overview of our simulation results. 
We find that sonic instabilities are a ubiquitous outcome in our simulations and excite acoustic 
modes that radiate angular momentum into both the star and the disk. In §[4j we show that the 
acoustic modes observed in simulations have a dispersion relation that is well-described by a few 



simple modifications to the dispersion relation of the supersonic vortex sheet (Miles 1958; Gerwin 
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19681 |Belyaev fc Rafihov|2012| | . In |3} we show that radiation from k z = acoustic modes entirely 



accounts for the angular momentum transport in our simulations and discuss the implications of this 
result for semi- analytic models of the BL. Finally, in ^6] we discuss the effects of disk stratification in 
the z-direction on acoustic modes and angular momentum transport, and we find that stratification 
makes no fundamental difference to our results. 



2. Numerical Model 
2.1. Governing Equations 



For our simulations we use the Godunov code Athena (Stone et al 2008) to solve the Euler 
equations of fluid dynamics in cylindrical geometry with an isothermal equation of state. These 
equations are 

^ + V-(pv) = 0, (1) 

+ V • (pvv) + VP + = 0, (2) 

P = ps 2 . (3) 

Here, v is the velocity, p is either the density (3D) or the surface density (2D), P is the pressure, 
<E> is the gravitational potential, and s is the sound speed. We also use w throughout this work to 
denote the cylindrical radius. 



2.2. Nondimensional Units 

Our initial setup, described in more detail below, consists of a star which attaches to a constant 



density disk via a thin interfacial region. Following the convention of Belyaev et al. (2012), we 
nondimensionalize quantities by setting the radius of the star to vo* — 1, which defines a unit of 
length, the Keplerian angular velocity to VLk{vo*) — 1, which defines a unit of time, and p = 1 within 
the disk, which defines a characteristic density. Expressed in these units, the orbital timescale at 
the surface of the star is Pk — 2tt. 

We also define a characteristic Mach number as M = 1/s, which is the true Mach number 
of the flow in the disk at the surface of the star, since Vk(^) = 1- For a given set of initial 
conditions, the Mach number is the sole physical parameter that determines the time-evolution of 
the system ( [Belyaev' etHp012| |, assuming an isothermal equation of state. Thus, the solutions to 



equations ([TJ-Q form a one-parameter family in the Mach number. However, in practice there are 
also numerical parameters, e.g. the size of the simulation domain, that influence the outcome of a 
simulation. 
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2.3. Gravitational Potential 

We take the gravitational potential to be fixed throughout the course of a simulation. This 



amounts to ignoring the self-gravity of the outer layers of the star (Cowling's approximation, Cowl- 
ing ( 1941| )), which was discussed in the context of the BL in Belyaev & Rafikov (2012). 



We run both vertically stratified and unstratifled simulations. For vertically unstratified runs, 
we use the cylindrically symmetric potential 



1 

w 



whereas for vertically stratified runs, we use the point mass potential 

* = - 1 



Vw 2 + z 2 ' 



(4) 



(5) 



where z is the height above the midplane. 



2.4. Initial Rotation Profile 

Initially, the star attaches to the disk via a thin interface, inside of which the rotational velocity, 
va, rises very nearly linearly from zero (the velocity in the star) to the Keplerian orbital velocity 
(the velocity in the disk). The initial rotation profile is given by 





ft(w) = < w' 



-3/2 
-3/2 



VJ—1 I 1 
S B L,0 ^ 2 



w<l- ™ "star 
1- S -^<uj<1+ S -^ 



'interface" 



(6) 



w > 1 + 



5bl,o 
2 



"disk" 



and is the same as in Belyaev et al. (2012) (see their Figure 1). For all our runs, we use 5bl,o = -01, 
so the interface is as narrow as possible, while still being resolved with ~ 10 cells. Our results 
are independent of 5bl,o, since over the course of a simulation, the interface thickens into a self- 
consistent BL with width ^> 5bl,o due to angular momentum transport by acoustic modes. 



2.5. Initial Density and Pressure Profiles 



The initial pressure profile is specified everywhere throughout the domain through the equation 
of hydrostatic equilibrium 



1 dP 

p dm 



dw 
ldP 

p dz 



dz ' 



(7) 
(8) 
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The density is directly proportional to the pressure and is related to it through the equation of 
state 

For the 2D and 3D unstratifled simulations, the initial density is p{w^ t = 0) = 1 within the disk, 
and p(zu, t = 0) = exp[— &(w)/s 2 } within the star, where O = 0. For the stratified simulations, these 
expressions are also accurate for the midplane density, po(zu, t = 0) = p{vo, z = 0, t = 0), as long as 
we use the midplane potential <&o(m) = ${w,z — 0). The initial density throughout the simulation 
domain for a stratified simulation can be expressed in terms of the initial midplane density as 
p(vj,z,t = 0) = po(vu,t = 0) exp[— ($(n7, z) — <f>o(w))/s 2 ]. This formula is valid everywhere, 
including the interfacial region, since we have rotation on cylinders throughout the simulation 
domain, i.e. ft is a function only of w (equation (6l). 



2.6. Boundary Conditions 

In all of our simulations, we use periodic boundary conditions (BCs) in the (f) direction and 
"do-nothing" BCs in the w direction. The "do-nothing" BC simply means that all fluid quantities 
retain their initial, equilibrium values for the duration of the simulation at the inner and outer 
cu-boundaries. A "do-nothing" BC in w is advantageous to an outflow BC, especially at the inner 
edge of the domain inside the star, since it prevents a systematic change in the density profile 
over the course of the simulation due to small deviations from equilibrium. A "do-nothing" BC is 
also advantageous to a reflecting BC, since it allows sound waves to leave the domain. However, 
a "do-nothing" BC does not absorb radiation perfectly (neither does an outflow BC), and there 
is some partial reflection of waves from the radial boundaries, although the reflection coefficient is 



significantly less than unity (see £5.1.2) 



For the 3D unstratifled and stratified simulations we use periodic BCs in the z-direction. The 
reason for using periodic BCs in the z-direction as opposed to outflow BCs, even in the stratified 
case, is that it is hard to achieve a good numerical equilibrium with outflow BCs. The application 
of outflow BCs is potentially more straightforward in spherical geometry (in the ^-direction), since 
in that case the gravity vector is parallel to the ^-boundary. 



2.7. Initial Perturbations and Burn-In 



Although we initialize our density and pressure profiles to satisfy the equations of hydrostatic 
equilibrium (equations [7] and [§]), we find that we need to burn-in our simulations for t ~ 100 
to let any initial transients damp out and achieve a good numerical hydrostatic equilibrium. This 



burn-in procedure was previously used by Armitage (2002) in studying the BL, and we find it easier 



than setting up an exact numerical equilibrium ( Zingale et al.|[2002 ). 

After the burn-in period, we introduce random perturbations to v w of amplitude Sv v 



.001 
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in the region w > 1. These random perturbations to the radial velocity seed the sonic instability 
and facilitate the growth of acoustic modes. The reason for only introducing perturbations in the 
region w > 1 is that waves traveling from high density to low density regions amplify by a factor 
of p -1 / 2 to conserve energy. Since the density in the star is up to 10 6 times higher than that in the 
disk, we want to avoid any transient waves that amplify as they travel out of the star and into the 
disk. 



2.8. Simulation-Specific Parameters 

Table [l] lists the simulation specific-parameters for each of our runs. Typically, our simulations 
are run until a time of t w 400 — 600 (70 — 100 orbits at the inner edge of the disk), where t — 
corresponds to the time when the random perturbations are introduced to the radial velocity. We 
find that this is long enough for the boundary layer to widen to an approximately steady-state 
width, after an initial period of rapid growth. Some of our simulations are run for over 100 orbits, 
and we observe little change in the behavior of the fluid in the vicinity of the BL or in the width 
of the boundary layer at the end of these runs as compared to our typical ones. 

The four most important parameters we vary from simulation to simulation are dimensionality 
(2D or 3D), Mach number (M = 6 or M — 9), stratification in the z-direction (unstratifled or 
stratified), and the <j) extent of the simulation domain. Comparing the effects of dimensionality 
and stratification across simulations with all other parameters held constant allow us to probe how 
the thickness of the disk affects the dynamics of the system. Varying the Mach number gives us 
a handle on how the phenomena we discuss translate over to astrophysical systems which tend to 
have higher Mach numbers (e.g. M ~ 100 for white dwarfs). Finally, changing the <j) extent of 
the simulation domain gives us a measure of control over the wavelength of the dominant modes 
present in the box, allowing us to probe the dispersion relation for our model (Q. 

3. General Overview 

We now present a general overview of our simulation results that will frame the discussion 
throughout the paper. 

3.1. Notation for Measured Quantities 

We begin by introducing some notation, regarding how simulation quantities are measured. 
For 3D simulations, we define a z-averaged density as 



E(zx7,<£) = 
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label 


dimension 


M 


stratified 


tu-range 


<^-range 


z-range 


N m xNjX N z 


3D6a 


3D 


6 


no 


(.7, 2.5) 


(0, 2tt) 


(-1/12,1/12) 

\ / > / / 


2048 x 2048 x 32 


3D6b 


3D 


6 


no 


(.7, 2.5) 


(0, 2tt) 


(-1/6,1/6) 


2048 x 2048 x 64 


3D6c 


3D 


6 


yes 


(.7, 2.5) 


(0, 2tt) 


(-1/3,1/3) 


1536 x 1536 x 128 


3D6d 


3D 


6 


no 


(.7, 2.5) 


(0, tt/3) 


(-1/12,1/12) 


2048 x 384 x 32 


3D9a 


3D 


9 


no 


(.85, 2.5) 


(0, 2tt) 


(-1/18,1/18) 


2048 x 2048 x 32 


3D9b 


3D 


9 


no 


(.85, 2.5) 


(0, 2tt) 


(-1/9,1/9) 


2048 x 2048 x 64 


3D9c 


3D 


9 


yes 


(.85, 2.5) 


(0, 2tt) 


(-2/9, 2/9) 


1536 x 1536 x 128 


3D9d 


3D 


9 


no 


(.85, 2.5) 


(0, it/ 3) 


(-1/18,1/18) 


2048 x 384 x 32 


3D9e 


3D 


9 


no 


(.85, 2.5) 


(0, 2tt/7) 


(-1/18,1/18) 


2048 x 384 x 32 


3D9f 


3D 


9 


no 


(.85, 2.5) 


(0, 2vr/7) 


(-1/18,1/18) 


2048 x 384 x 64 


2D6a 


2D 


6 




(-7, 2.5) 


(0, 2vr) 




2048 x 2048 x 1 


2D6b 


2D 


6 




(-7, 2.5) 


(0, tt/2) 




2048 x 512 x 1 


2D6c 


2D 


6 




(-7, 2.5) 


(0, vr/4) 




2048 x 256 x 1 


2D6d 


2D 


6 




(-7, 2.5) 


(-.4, .4) 




2048 x 512 x 1 


2D9a 


2D 


9 




(.85, 2.5) 


(0, 2tt) 




4096 x 4096 x 1 


2D9b 


2D 


9 




(.85, 2.5) 


(0, tt/2) 




4096 x 1024 x 1 


2D9c 


2D 


9 




(.85, 2.5) 


(0, tt/4) 




4096 x 512 x 1 



Table 1: Summary of simulation parameters. The columns from left to right are: simulation label, 
dimensionality (2D or 3D), Mach number, stratification (yes = stratified, no = unstratifled), radial 
extent of the simulation domain, azimuthal extent, vertical extent (only in 3D), number of grid 
points in the radial x azimuthal x vertical directions. 



where the integral is performed over the entire z-extent of the domain, and Az is the thickness of 
the simulation domain in the z direction. The ^-averaged density and perturbations to the density 
are then defined as 

EoM = ^ y d0E (10) 

5S(tU,0) = E-E (11) 
5p(w,(f),z) = p-po- (12) 

Next, we define the density- weighted averaging operator represented by a set of triangle brackets 



as 



'^A^/"^ (13) 

In 2D, the density- weighted average is the same, except that there is no integral over the z- 
dimension. Sometimes, we perform density-weighted averages in only the z-direction so we also 
define the operation 

(/> = = ArE /***• <"> 
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Density- weighted averages are the proper way of measuring simulation quantities, since math- 



ematically exact analogs to the fluid equations can be formulated in terms of them (Balbus & 



Papaloizou 1999). For instance, one-dimensional forms for the continuity and angular momentum 



transport equations are 



aSo , 1 9 ( v / U n 
-^-^ — (iuSo(^)) 

at w aw 
d 1 d 

-( zu 2 z n) + -—(w 2 x (v 4> v^)) = o, 



(15) 
(16) 



where we have assumed molecular viscosity is negligible in equation (16), and we have defined the 
angular velocity as 



Q = 



(17) 



Equation (16) expresses angular momentum conservation and the second term is the divergence 



of the angular momentum flux, 

F L = w^ (v^)- (18) 

Throughout the paper we shall also refer to the angular momentum current, which is denned as 

C L = 2tvzuF l . (19) 

The angular momentum current has the intuitive interpretation that the rate of change of angular 
momentum in an annular region is the difference in angular momentum current between the two 
edges of the annulus. 



The angular momentum flux, Fl, can be decomposed into two parts (Balbus & Papaloizou 



1999), which are the advective and stress terms 

F L = F A + F S 
F A = w 2 ^{v w ) 
F s = w^ (5v^v m ) 
Sv^ = — wft. 

These definitions are trivially generalized to 3D case by replacing with po- 



(20) 
(21) 
(22) 
(23) 



The advective term corresponds to angular momentum transport by the bulk motion of the 



fluid and can be written in terms of the mass accretion rate, M, as 



F A = 



M 



-ZD 



2 Q 



2nvu 

M = -2ttzd^q(v w ) 



(24) 
(25) 
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The stress term, on the other hand, corresponds to angular momentum transport by turbulence 
and/or waves. In the standard o:-viscosity treatment, the stress term is alternately written as 
( |Shakura fc Sunyaev 1 |1973| |Popham fc Narayan|[l995l |Balbus fc Papaloizou|[l999l ) 



F S 
F S 

F S 



2 ( dlnfi 
\ dmw 

2 / dlnfi 

map^os -— 



1 + 



s d\nP 
ft dzu 



-1/2 



(26) 
(27) 

(28) 



For any of these (or any other) a viscosity model to provide a good description of the BL, a should 
be approximately independent of w. However, all conventional a viscosity models are based on the 
notion of local turbulent stresses. As we discuss in in our BL simulations angular momentum 
transport is facilitated by waves not turbulence, and so is not well described by any local model. 



3.2. Simulation Results 



We now paint a picture of the typical evolution of one of our simulations. Due to the large 
shear in the interfacial region (j|2.4|), our initial rotation profile is unstable to the sonic instability 



(Glatzel 1988; Belyaev &; Rafikov 2012; Belyaev et al. 2012[ ). Consequently, acoustic modes are 
excited starting from the initial seed perturbations to the radial velocity (j |2.7| ). The amplitude of 
these acoustic modes increases with time until t <~ 60 when the sonic instability saturates. 

Around t ~ 60, there is typically a single mode present with a well-defined azimuthal pattern 
number, m, and pattern speed, ftp. This mode, which is discussed in detail in §4.2. 1\ has ftp ~ 1 — 5 



and radiates angular momentum in the form of sound waves away from the vicinity of the interfacial 
region into both the disk and the star. This radiation of angular momentum causes the interfacial 
region to widen and drives accretion of material from the innermost portion of the disk onto the 
star. 

Fig. [T^i shows a radius-time image of for simulation 3D9a. From the image, it is clear that 
a gap in the density is opened up in the innermost portion of the disk starting at around t <~ 60, 
which is around the same time when the sonic instability saturates. The resulting radial density 
gradient gives rise to partial pressure support for fluid in the inner disk and causes the rotation 
profile to deviate from Keplerian, an issue that we return to in §5.4| Figs. [2k & b show plots 



of and mfl, respectively, at various times for simulation 3D9a. It is evident that as the inner 
region of the disk is depleted of mass, the velocity profile becomes increasingly sub-Keplerian. Both 
effects can be attributed to radiation of angular momentum away from the inner disk which drives 
accretion onto the star. 



Fig. [TJd shows a radius-time image of the radial velocity, and it is clear that there is radial 
infall in the inner region of the disk, which occurs during the same time as depletion of the inner 




Fig. 1. — Panels a Sz b show radius-time plots of and (v m ), respectively for simulation 3D9a. 



disk in Fig. [IJi. The striations in Fig. [TJ) correspond to density waves that emanate away from 
the boundary layer, and their "slant" indicates that by causality they are indeed emitted from the 
vicinity of the BL. 

It is evident from Fig. that after a period of vigorous infall, the radial velocity and hence 
the accretion rate in the inner disk drops off around t ~ 250. This can be attributed to a reduction 
in the amplitude of waves emitted from the BL, so the rate of angular momentum transport is 
diminished. Moreover, the pattern speed drops significantly during this period, and it is possible 



for trapped shocks to develop between the BL and an evanescent region in the disk (Belyaev et 



al. 2012). The reason for this shift in behavior seems to be related to the angular momentum 
budget of the inner disk, and we postpone a more detailed discussion of this issue until where 
we investigate angular momentum transport in depth. 



Finally, we comment on a result of Belyaev et al. (2012), where they found the width of the 
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Fig. 2. — Panels a & b show Sq and wVt vs. time, respectively, for simulation 3D9a. 



BL to extend across 6 — 7 scale heights in the 2D hydro case. Defining the BL as the region in 
which 

0.1 < (v^(w))/v K (w) < 0.9, (29) 

we find general agreement with the 2D results of |Belyaev et al.| Q2012| ) in our 3D simulations. 

Fig. [3] shows the width of the BL (in scale heights (equation [36])) as a function of time for 
both our 3D unstratifled simulations that span the full 2tv in azimuthal angle (solid lines) and the 
3D stratified simulations (dashed lines). The color denotes Mach number with black corresponding 
toM = 9 and red corresponding to M = 6. 

We point out that there is a sizeable degree of scatter among the simulations in Fig. [3j but 
that the mean is around 6 scale heights, which is consistent with 2D simulations. One possible 
reason for the larger amount of scatter in the 3D case is that the 2D simulations were run for a 
longer time, and the measurement of BL width was performed at a later time in the runs. Thus, 
the BL in our 3D simulations may still be settling down to its steady state width. 

We now turn to the modal structure of the BL, since it is ultimately waves which transport 
angular momentum in our simulations. Thus, understanding the dispersion relation of the BL is 
instrumental to understanding angular momentum transport by these waves. 



4. Boundary Layer Modes 



We find in our simulations that the shear in the boundary layer can generate three distinct 
types of k z = acoustic modes. These three types of modes are directly related to the upper, 



10 




time 

Fig. 3. — Width of the BL in scale heights as a function of time. The dashed curves correspond to 
3D stratified simulations and the solid lines to 3D unstratifled simulations spanning the full 2tv in 
azimuthal angle. The color differentiates Mach number with M = 9 in black and M = 6 in red. 



lower, and middle branches discussed in Belyaev & Rafikov (2012) for a 2D plane parallel vortex 
sheet in the supersonic regime. However, the dispersion relation for these modes is modified in a 
simple way by the radial stratification in the star and the Coriolis force in the disk. To begin our 
discussion, we first review the dispersion relation of the 2D plane parallel vortex sheet. 



4.1. Plane Parallel Vortex Sheet 

In the absence of stratification or Coriolis force, the dispersion relation for a 2D vortex sheet 
(with x and y coordinates locally representing w and <j) coordinates in the BL) with velocity profile 



sound speed s + for x > 0, and sound speed S- for x < in the limit V ^> s is given by (Belyaev 
fc Rankov||20l2l ) 

{V — s+, Upper branch 

V (j^r) , Middle branch (31) 
«s_, Lower branch. 

Note that by pressure balance p+/p~ = (s_/s + ) 2 , assuming 7 is everywhere constant. All three 



branches in equation (31) are sonic in nature, with outgoing sound waves emitted from the vortex 
sheet at x = 0. 



The upper and lower branches have the property that the sound waves emitted from the vortex 
sheet in the half-planes x > and x < 0, respectively, propagate nearly parallel to the vortex sheet. 
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In other words, \k y /k x \ ^> 1 for the upper and lower branches in the half-planes x > and x < 
respectively. The middle branch has the property that \k y /k x \ is the same on both sides of the 
vortex sheet. Thus, the wavefronts form the same angle with the vortex sheet in both the x > and 
x < half-planes. This is true even if the densities are different in the two half-planes. These sets 
of properties for the isothermal vortex sheet are instrumental in determining how the dispersion 
relations for each of the three wave branches are modified in going to the disk-star system. 

Panels a-c of Fig. [4] are schematics, which depict the geometry of the wavefronts for 8v Xl the 
perturbed velocity perpendicular to the vortex sheet, for each of the three branches. The arrows 
show the direction of the flow, and the solid vertical line is the vortex sheet at x = 0. At the 
vortex sheet, 5v x undergoes both a change in amplitude and a change in sign (see equation (36) of 



Belyaev & Rafikov (2012)). Thus, the wavefronts of 5v x appear shifted in phase by tt across the 
vortex sheet. 
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Fig. 4. — a) Wavefronts of 8v x for the lower wave for M = 6. The arrow depicts the direction of 
the flow in the x > half-plane. The vortex sheet itself is depicted by the solid vertical line and is 
located at x = 0. The solid lines depict wave crests, and the dashed lines depict troughs. Notice 
that the perturbation to Sv x undergoes a change in sign (phase shift by tt) at the vortex sheet. 
Panels b) and c) are the same as a), but for the middle and upper waves respectively. 
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4.2. Dispersion Relations for the Isothermal Boundary Layer 

We now discuss how the dispersion relations for each of the wave branches in the isothermal 
vortex sheet given by equation (31) are modified in going to the disk-star system. As a result, we 
find accurate, though, heuristic dispersion relations for the generalization of these modes to the 
isothermal boundary layer studied in our simulations. 



4-2.1. Upper Branch 



As discussed in £4.1, the upper branch for the vortex sheet is characterized by a sound wave 
with phase speed V — s+ in the half-plane x > 0, which propagates nearly parallel to the vortex 
sheet, i.e. \k y /k x \ ^> 1. 

In order to generalize this to the disk-star system, we postulate that in the disk-star system 
there is a branch characterized by a density wave in the disk that has {k^/k^l ^> 1 at w — w*. 
Here, and throughout the paper, we shall define k§ = m/vj, where m is the azimuthal pattern 
number. 

Since, {k^/k^l ^> 1 in the disk near the BL, the wave fronts are very loosely wound, and we 



cannot apply the classical dispersion relation for a fluid disk in the tight- winding limit (Goldreich 



fc Tremaine||1978 ). However, it is straightforward to derive an approximate dispersion relation for 
a loosely wound density wave, if we additionally assume m^$> 1 (Appendix [A]). The result is 



(ft - ft P ) 2 



+<-) ■ 



(32) 



Equation (32) generally has two solutions with one having ftp > ft and the other ftp < ft. For 



the isothermal vortex sheet, it turned out that the physical solutions, corresponding to outgoing 
waves have a phase velocity uo/k < V ( Belyaev fc Rafikov|2012[ ). By analogy, we take the solution to 
equation (32) having ftp < ft to be the physical one. In that case, taking w — zu+, V = m+ftfa*), 
and making the associations uj — mftp, and k y = m/m*, we see that equation (32) is analogous to 



equation (31) for the upper branch of the vortex sheet, except for the additional term containing 



the epicyclic frequency, which is due to the Coriolis force. 

We can generalize equation ([32]) to arbitrary values of k w by writing 



m 2 (fl - ftp) 2 = s 2 (kl + kl) + k 2 



(33) 



This approximate dispersion relation reduces to the one for tightly-wound waves in the limit 
Ik-cu/k^] ^> 1 (Goldreich & Tremaine 1978) and to equation (32) in the limit Ik^/k^] ^> 1. Thus, 
it is accurate in both the tight-winding and loose-winding limits. The essential modification to 
the tightly-wound dispersion relation is the replacement of the radial wavenumber with the total 



wavenumber inside the braces in equation (33) i.e. k. 



,2 
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An important feature of equation (33) is that it predicts an evanescent region in the disk 



directly adjacent to the BL. This can be deduced by noting that the left hand side of equation (33) 
vanishes when £1 = ftp, implying that k m is imaginary at the corotation radius. 

Table [2] compares the value of ftp observed in simulations when the upper wave is present to 
that predicted by the dispersion relation (32). For simplicity, we set w — and use the Keplerian 
values of ft and k when calculating the predicted values from equation (32), i.e. all variables 



are evaluated just outside the BL. Nevertheless, the agreement between the analytical predictions 
and the simulations is good to the level of several percent for both 2D and 3D simulations and 
irrespective of stratification. However, some of our theoretical estimates tend to be slightly high, 
especially at late times {t > 60 orbits) when the boundary layer has thickened substantially. One 
possible cause for the slight deviation between the measured and analytically predicted pattern 
speeds is that the boundary layer is not inflnitesimally thin, but we have assumed so in evaluating 
the dispersion relation at w 



label 


time 


m 


ftp measured 


ft F 


predicted 


2D6a 


20 


« 40 


.84 




.831 


2D6b 


20 


« 40 


.835 




.831 


2D6c 


20 


« 40 


.835 




.831 


3D6c 


120 


9 


.77 




.80 


3D9a 


100 


10 


.83 




.850 


3D9d 


60 


12 


.85 




.861 


2D9a 


30 


« 29 


.87 




.884 


2D9a 


280 


11 


.815 




.856 


2D9b 


30 


« 24 


.87 




.881 


2D9c 


30 


^> 32 


.88 




.885 



Table 2: Comparison between the analytical value of ftp and the value measured in simulations for 
the upper branch. The columns are from left to right: simulation label, time since the start of the 
simulation, m, pattern speed measured from the simulation, pattern speed predicted from equation 



(32). 



Table [2] shows that the behavior of the phase speed of the upper mode is in reasonable agree- 
ment with the simple vortex sheet prediction (31). Indeed, the latter would result in ftp = 
uojikyW*) « 0.83 and 0.89 for M — 6 and 9 correspondingly. This agrees with the trend seen 
in our simulations, where a higher value of M results in a somewhat larger value of ftp. Also, the 



magnitude of ftp is close to f2x(^*)> m accord with equation (31) for high Mach number. 
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4-2.2. Lower Branch 



For the case of the isothermal vortex sheet (j |4.1| ), the lower branch is characterized by a sound 
wave which propagates nearly parallel to the vortex sheet in the half-plane x < with phase velocity 
S-. We now discuss the modifications that need to be made to this dispersion relation in order to 
apply it to the disk-star system. 

We begin by postulating that the lower branch in the disk-star system is characterized by a 
sound wave which propagates azimuthally in the star. One immediate difference as compared to 
the isothermal vortex sheet case is that the star is stratified, so the phase velocity of the lower 
branch is not simply s, but is modified by the effects of gravity and stratification. 

In a non-moving, plane-parallel atmosphere stratified in the x-direction with constant sound 
speed, «s, and gravity, g, the dispersion relation for gravito-sonic waves (p and g modes) is given by 



e.g. |Vallis | (|2006j) as 



UJ 



2 



1± 1 



4(7~l)fcg 



I 2 

h s 



k 2 x + k 2 y + l/Ahl 



(34) 

(35) 
(36) 



where h s is the pressure scale height. Note that the dispersion relation (34) applies to all wave 



lengths, not just to those that are short compared to the scale height of the atmosphere. 



The upper sign inside the braces on the right hand side of equation ([34]) corresponds to p- 
modes and the lower sign to g-modes. An isothermal equation of state (7 = 1) is neutrally buoyant, 
so g-modes have uj = 0, and we will not consider them further. On the other hand, the dispersion 
relation for p-modes with 7=1 reduces to 



UJ 



1 



s~ I k x + k y + 



(37) 



Note that the effect of gravity and stratification enters into the dispersion relation through the 
1/4/ij? term, and in the limit h s 00 we recover the dispersion relation for a sound wave in a 
uniform medium. However, since we typically find A > h s in our simulations for the dominant 



modes, the stratification term in equation (37) is very important for comparing with simulations. 



We can use equation (37) to derive a dispersion relation for the generalization of the lower 



branch to the disk-star system. We start by setting k y = m/m, and by analogy with the vortex 
sheet, we assume that \k y /k x \ ^> 1 inside the star. Making again the association, uj = mf2p, we 
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have 



n P 




+ 



1 



2mh s (wo) 



= Q(ro*) — < 



2m vjq 



(38) 
(39) 



Since ftp must be independent of radius, we have chosen an "effective radius", zuq < zu* 3 to use in 



the expression on the right-hand side of equation (38). Although the choice of an effective radius 



may seem like an "ad-hoc" prescription, we provide theoretical justification for it in Appendix [Bj 
It is further justified by the agreement between the value of fip(m, M) measured in simulations 



and that given by equation (38). 



Figure [5] shows a plot of Op vs. rn from our simulations when the lower wave is observed. It is 



evident that the agreement between equation (38) and the data is quite good over a wide range of 



values for rn and ftp. Moreover, the 2D simulations (circles), 3D unstratifled simulations (squares), 
and 3D stratified simulations (triangles) lie on the same curve for both M = 6 and M = 9, implying 
that the behavior is fundamentally the same in both 2D and 3D and with or without stratification. 

It is of note that the lower mode in the BL has phase speed considerably higher than ftp = 



M 1 ft(m+) predicted by the simple vortex sheet model, see equation (31). This is because the 
results presented in §4.1| were obtained assuming an unstratifled star, in which case h s —> oc and 



our BL dispersion relation (38) would reduce to its vortex sheet analog. However, in an isothermally 



stratified star, the density gradient strongly modifies the dispersion relation, as evidenced by the 



second term in equation (39), and makes ftp not only a function of M but also of m, the azimuthal 
wavenumber. 



4.2.3. Middle Branch 

For the isothermal vortex sheet with a constant sound speed everywhere throughout the domain 
(s+ = the phase velocity of the middle branch is simply V/2 (j |4.1| ). Moreover, the wavefronts 
make the same angle with respect to the vortex sheet both above and below it, which is true even 
for the case of non-equal densities above and below the vortex sheet (|Belyaev fc Rafikov 2012). 



This last fact can be expressed mathematically as k Xj +/k y = k x -/k y , where the plus and minus 
signs denote the value of k x in the x > and x < half-planes, respectively. 

We now postulate that for the disk-star system the wavefronts in the star and in the disk form 



the same angle with respect to the azimuth across the boundary layer. Using equations (33) and 
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Fig. 5. — Comparison between the analytical dispersion relation, fip(m, M), and simulation results 



for the lower branch. The solid line is the analytical prediction from equation (38) for M = 6 using 



an effective radius of wo — .84, and the dashed line is for M = 9 using wo = .90. The points show 
values of Op and m measured from simulations when the lower branch is dominant. The circles 
represent 2D simulations, the squares 3D unstratified simulations, and the triangles 3D stratified 
simulations. Filled points represent simulations spanning the full 2tv in azimuthal angle, and open 
points represent simulations for a wedge which spans an azimuthal angle < 2tt. There is a solid 
triangle and a solid square that have merged near the solid curve at M = 6. 



(37), we can write 



zi7,disk 



^tu.star 



(u> - mn) 2 + k 2 
<M 2 



CO 



(sk^) 2 {2h s k <j) f 



(40) 
(41) 



Setting equations (40) and (41) equal to each other at w — (the condition that the wavefronts 
in the disk and star form the same angle with respect to the azimuth across the BL), the dispersion 
relation for the middle branch is 

o2 



1 + 



1 



4«* 



(42) 



Table [3] lists the pattern speeds and m-numbers measured in simulations when the middle 
branch is clearly dominant and compares the pattern speed to that predicted by equation (42). For 
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simplicity, we again assume a Keplerian profile for k and Vt when computing the predicted pattern 
speed. The agreement between the predicted and measured values is at the several percent level. 



Also, Vlp is close to f2(cc7*)/2, in agreement with the vortex sheet model, see equation (31). The 
discrepancy between the measured and analytically-predicted pattern speeds could arise because 
the condition k^^/k^ = k^-jk^ is not exactly satisfied for the disk-star system as it is for the 



isothermal vortex sheet. Nevertheless, we find that equation (42) provides a good approximation 
to the pattern speed observed in our simulations when the middle branch is dominant. We also 
note that we have fewer measurements for the middle branch as compared to the upper and lower 
branches, because the middle branch is a transient phenomenon observed only during the linear 
growth regime of the sonic instabilities, if at all. 



label 


time 


m 


ftp measured 


Qp predicted 


3D6a 


10 


« 42 


.52 


.502 


3D6d 


10 


« 36 


.515 


.503 


3D6c 


10 


« 44 


.52 


.502 


3D9e 


20 


r& 28 


.57 


.515 



Table 3: Comparison of the analytically predicted value of fip, and the value measured from sim- 
ulations for the middle branch. The columns of the table are from left to right: simulation label, 
time since the start of the simulation, m, pattern speed measured from the simulation, pattern 



speed predicted by equation (42). 



4.3. Spatial Morphology of the Upper, Middle, and Lower Branches 

We have identified three branches in our simulations that are the generalizations of the upper, 
middle, and lower branches of the 2D isothermal vortex sheet. We now discuss the structural 
morphology of each of these branches. 

Fig. ^ shows snapshots of vj\/Yj{v V j) z vertically-averaged in the z-direction for each of the 
three wave-branches from 3D unstratifled simulations. The quantity w\/T^{v vu ) z is approximately 
constant throughout the star and disk, so it is especially useful for visualizing the radiation pattern. 
This is because waves conserve energy and angular momentum current, which are both of order 
vo 2 Yi5v 2 , where 8v is the magnitude of the velocity perturbation. The dashed and dotted vertical 
lines in Fig. [6] depict the edges of evanescent regions and corotation radii respectively. Notice that 
there are always two corotation radii: one inside the boundary layer at vo w 1, and one inside 
the evanescent region within the disk. In principle, there is also an evanescent region around the 
corotation radius inside the BL, but this region is narrow, due to the steep radial gradient of angular 
velocity within the BL. 

The sharp features seen in this Figure (especially in panel b) are the shocks into which the 
propagating acoustic modes evolve as a result of their nonlinearity. Even though the nonlinearity 
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1 1.2 1.4 1.6 1.8 2 2.2 2.4 

radius 

(b) 




0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 

radius 



(c) 

Fig. 6. — Panels a, b, c show snapshots of w\/Yi{v vu ) z , when the upper, lower, and middle branches 
are dominant. Panel a corresponds to simulation 3D9d at t — 65, b to simulation 3D9e at t — 320, 
panel c to simulation 3D6a at t — 10. Dashed vertical lines depict the boundaries of evanescent 
regions and dotted vertical lines depict corotation radii. 



is relatively mild, with v w <~ 0.1s, implying density jump of ~ 10% across the shock fronts in this 
particular simulation (see also Fig. 10), the waves propagating into the disk break and become weak 
shocks. In our simulations we typically find this to occur not too far from the BL, at separations 
of (0.2 — \)vo* away from their launching point. 
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The solid and dashed curves in Figure [6] trace the shape of shock fronts or wave fronts if the 
wave has not shocked. The formula for the shape of a front is simply 

d$ = ± k 
dw ka. 



where k m is given by formula (33) inside the disk or (37) inside the star. The sign on the right 
hand side of equation (43) determines whether the wave front is angled up or down, assuming we 
take k m and k^ to be positive. We see from Figure [6] that the analytic curves for the shock fronts 
are in excellent agreement with snapshots from the simulations, which is due to the fact that waves 
are only weakly nonlinear. 

For the upper branch, which has high ftp close to unity, the evanescent region extends all the 
way down to the boundary layer, so density waves emitted from the BL that propagate into the 
disk have a pattern speed that is faster than the angular speed of the disk material and wind up 
with increasing radius. For the lower and middle branches, however, the evanescent region in the 
disk is well-separated from the boundary layer, so density waves emitted from the BL have a slower 
pattern speed than the angular speed of the disk material. These density waves unwind as they 



travel towards the forbidden region, partially reflecting off it, with very little transmission (see {5.2 
for a lower bound on the reflection coefficient measured from simulations). 

As a result of this reflection, a standing shock pattern can be set up in the disk between the 
BL and the forbidden region (Fig. ^jp) if a resonance criterion is satisfied. This criterion states that 
in order to achieve standing shocks, the total azimuthal angle subtended by a shock in traveling 
from the BL to the forbidden region and back must be an integer multiple of the angle subtended 
by one wavelength of the mode. Mathematically, this criterion is 

^ = 2A0 shock (44) 
m 

A0 s hock / dw—, (45) 

<W dw 



see 



Belyaev et al. (2012) for more details. Here, wq is the inner edge of the forbidden region in 



the disk, wbl ~ 1 5 and d(j)/dw is given by equation (43). The integer p > 1 gives the number of 
shock crossings for a shock traveling from the BL to the forbidden region or vice versa, counting 
the reflection point at the edge of the evanescent region as a shock crossing. For instance, solving 



equation (44) for the m — 14 mode depicted in Fig. |6p we find p = 6.03. Thus, p is very close to 
an integer, which allows the standing shock pattern to form, and it is evident that each outgoing 
and incoming shock undergoes six shock crossings (counting the point of reflection). 

Although p does not necessarily need to be an integer, it is likely that modes with an azimuthal 
wavenumber, m, that is close to an integer value of p are reinforced relative to modes with a similar 
value of m, but for which p is not close to an integer by the resonance condition. 
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4.4. Modes with k z ^ 

The upper, lower, and middle branch modes are all two-dimensional in the w — (j) plane in 
the sense that they do not require the z-dimension to operate. In principle, there are also Kelvin- 



Helmholtz (KH) modes for vertical wavenumbers kz/k^ > M (Belyaev & Rafikov 2012). 



To determine the relative importance of the k z ^ vs. the k z = modes, we compute the 
power spectrum from the unstratified simulations along the z-dimension. We work in the variable 



5p/v£o, which by the definition of 5p (equation 12 ) sums to zero when integrated over the entire 



simulation domain. This is a natural variable to use, since one expects the time-averaged value 
of \5p/\/T^o\ to be roughly constant for a wave across the entire simulation domain up to factors 
related to the geometry and rotation profile. 

To obtain the power spectrum, we first perform a discrete Fourier transform in the z-direction. 
In order to carry out summations, we let Sp/y/Eo be a function of the grid points (/,m,n) in the 
(m, 0, z\ directions respectively. The Fourier transform in the z-direction then becomes 

where N z is the number of cells in the z-dimension. The power spectrum is then defined as 

N^-l^-l 



P n = — — Yl J2^n(l,m)) 2 - (47) 



/=0 m=0 

We note that by Parseval's theorem, 

N z -1 . N z -lN^-l N^-l 



E'. = ^EEE ^ • («) 



Fig. [7| shows power spectra computed from the simulations. Plotted is log 10 (P n / ^ n P n ) vs. 
k Zl where k z = 27rn/Az, and Az is the width of the simulation domain in the z-direction. The 
black lines are from simulation 3D9e, and the red lines are from simulation 3D9f which has the 
same parameters but twice the z-resolution. The solid lines correspond to a time when the upper 
mode is dominant, and the dashed lines to a time when the lower mode is dominant. From the 
figure, it is clear that the n — mode has at least several orders of magnitude more power than all 
of the n > modes combined. Note that we only plot the modes from n — to n — N z /2, where 
N z /2 = 16 for simulation 3D9e. This is because modes having n — N z /2 + 1 to N z — 1 (if N z is 
divisible by 2) are just the complex conjugates of modes having n — N z /2 — 1 to n = 1, for a real 
signal with no imaginary component, such as ours. Simulation 3D9f has N z /2 = 32, but we only 
show the modes up to n — 16, since there are no interesting features past n > 16, and the power 
spectrum simply continues to decay with increasing n. 
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It is clear from Fig. [7] that there is more power in the k z ^ modes relative to the k z = 
mode for the higher resolution simulation as compared to the lower resolution simulation. Thus, 
it is possible that our simulations are not fully converged in the z-direction. Nevertheless, the fact 
that the power in the k z ^ modes is still several orders of magnitude less than the power in the 
k z = mode even for the higher resolution simulation suggests that the k z = mode is indeed 
dominant irrespective of the resolution in z. 




_gl 1 1 1 1 1 1 1 1 1 1 

100 200 300 400 500 600 700 800 900 



Fig. 7. — The power spectrum, log 10 (P n / J2 n ^n) as a function of k z . The solid lines are for the 
upper mode and the dashed lines for the lower mode. Simulation 3D9e is shown in black and 3D9f 
in red. Note that simulations 3D9e and 3D9f have the same vertical extent, A z , but a different 
resolution in the z-direction. When comparing across simulations with different resolutions, one 
should compare k z rather than n, since modes with the same value of k z correspond to the same 
vertical wavelength, which is not true for modes with the same n if the simulations have different 
resolutions. 

Fig. [8]shows a fc^-time plot of log 10 (P n ) for simulations 3D9e and 3D9f. One can see interesting 
temporal features, such as rapid excitation of both k z ^ modes and k z — modes and their gradual 
subsequent decay. One possible source for these periodic outbursts are secondary KH instabilities 



Belyaev et al. (2012). However, in the 3D case, unlike the 2D case, a kink-like inflection point in 
the azimuthal velocity profile does not form, so this mechanism is not certain. It is also evident in 
Fig. [8] that the higher resolution simulation (3D9f) has a greater proportion of power in the k z ^ 
modes, but again the k z — mode is clearly dominant for the duration of the simulation even in 
the higher resolution case. Thus, one would expect the k z — modes to dominate the angular 
momentum transport in the system, which we show explicitly in the next section. 
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Fig. 8. — a) fc^-time plot of log 10 (P n ) for simulation 3D9e. b) k z -time plot of log 10 (P n ) for simulation 
3D9f, which has all the same parameters as simulation 3D9e, except that the z-resolution has been 
doubled. We again point out that comparing across simulations with different resolutions, one 
should compare k z rather than n, since modes with the same value of k z correspond to the same 
vertical wavelength. 



5. Angular Momentum Transport 

Here we discuss angular momentum transport due to the upper and lower wave branches. We 
find good agreement between theory and simulation and conclude that the upper branch is more 
efficient at transporting angular momentum. We do not consider the middle branch here, since it 
is dominant only for short periods of time while the instabilities are still growing. 

As we shall soon demonstrate, angular momentum transport in our simulations is dominated 
not by turbulent stresses, but rather by acoustic radiation of angular momentum away from the 
boundary layer both into the star and into the disk. Moreover, the k z = modes dominate the 
angular momentum transport, as hinted at in §4.4[ which allows us to use the theory of density 



waves in thin disks to calculate the rate of angular momentum transport (e.g. Appendix J of Binney 



fc Tremaine| ([20081)) 



The angular momentum current through radius w is given by equation (19). By analogy with 



equation (20), it can be decomposed into a stress term and an advective term. 



C s (vj) = 2ttvjF s 



(49) 
(50) 



Waves transport angular momentum exclusively through stresses, and we can use equations 



(49) and (22) to calculate Cs- Inside the star, Cs for a given mode with azimuthal wavelength m 



is given in terms of the amplitude of the density perturbation 5H m (iu), as 



rn 



n 2 P 



So 



(51) 



Here we have used the relations in Appendix fBl to express Sv^^ and Sv^^ in terms of 5S r 



^m- The 

total angular momentum current due to waves is the sum over all modes for both incoming and 
outgoing waves and is given by 



(52) 



We next assume that all waves propagate outwards from the BL, and that k m va/m is relatively 
insensitive to m near the dominant m value, which we denote as rn. This pair of assumptions will 



be verified later by comparison with simulation results. However, it fixes the sign in equation (51) 



and allows us to use Parseval's theorem ( Belyaev et al.|[2012 ) to write the total stress due to waves 
in the star as 

= - (him) glT «(¥)'■ <53) 



rn 



n 2 p 



The analog to equation (53) for the angular momentum current due to density waves in the 



disk is ( |Binney fc Tremaine||2008l |Belyaev et alT1[20T2| ) 



m 



(n/m) 2 - (ft - ftp) 



(54) 



where we have again used Parseval's theorem and assumed that all waves propagate away from the 
BL. 

Waves propagating into the disk shock, resulting in dCs/dzu ^ due to dissipation of the 
wave with increasing distance away from the BL. The shock-related change in angular momentum 
current carried by the wave as a function of radius in the disk is given by ( Larson|[l990 ; Belyaev et 
aLlf20T2l) 



dC s 
dw 



-mzuY,r 



dE 
' dm 



(55) 
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Here, dE/dm is the energy converted into heat per unit mass of fluid crossing the shock, which 
comes at the expense of the energy contained in the shock and causes it to damp. For the isothermal 



case, the expression of Larson (1990) for dE/dm is 

dE s 2 e 3 



dm 6 (1 + e) 2 



(56) 



in the limit e«l, where e = AS/S is the jump in density across the shock divided by the preshock 
density. 

It is also possible to derive an expression for dE/dm which is valid for arbitrary values of e in 
the isothermal limit 7 — )> 1 (Appendix [C]) 

dE_ 2 e(2 + e)-2(l + e)ln(l + e) 

dm~ S 2(1 + e) ■ {b<} 

This converges to Larson's result (equation [56]) to leading (third) order in e. 

In contrast to waves propagating into the disk, waves propagating away from the BL and into 
the star do not shock, because their amplitude is proportional to p -1 / 2 by conservation of energy 
and angular momentum. Thus, as waves propagate into the star in our simulations, their amplitude 
is reduced, nonlinear effects are mitigated, and dCs/dw — 0. Inside a real star, this reduction in 
amplitude could not continue indefinitely, and the waves would eventually be absorbed, depositing 
their energy and angular momentum to the stellar fluid. However this absorption could potentially 
occur deep inside the star, making the BL appear cooler than if the energy were dissipated locally. 

We now apply these results to understand angular momentum transport by the upper and 
lower acoustic modes, which appear in our simulations most often. 



5.1. Upper Branch 

5.1.1. Disk 

Fig. [9] shows Cs inside the disk for simulation 3D9d at t ~ 75, which corresponds to approx- 
imately the same time as Fig. [6k, when the upper mode is clearly dominant. The solid black 



curve shows the measured value of C#, using equations (22) and (49), and the red curves are the 



analytically-expected values of Cs, calculated using equation (54) together with the value of AE/E 
across shock fronts measured from simulations. The solid red curve uses the standard WKB ex- 
pression for k m (Goldreich & Tremaine 19781), and the dashed red curve uses k m from equation 



(33). The standard WKB expression can be obtained from equation (33) by eliminating the k, 



term in the brackets. For simplicity, we use the Keplerian values of n and in equation (54) when 



calculating the analytically-expected value of Cs (red curves). This is a good approximation, since 
the rotation profile is very nearly Keplerian in the disk, away from the direct vicinity of the BL. 
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Fig. 9. — Plot of Cs for simulation 3D9d at time t = 75 within the disk. See text for details. 



Regardless of which form for k w is employed, equation (54) itself is only valid in the tight- 
winding approximation. Thus, the red curves converge to the black at large but deviate signif- 
icantly as we approach the evanescent region around corotation, where the upper mode is loosely 
wound. For the solid red curve, the evanescent region is located between the vertical dashed lines, 
which correspond to the two Lindblad radii, and for the dashed red curve, the evanescent region is 
located between the dotted vertical lines. The corotation radius is marked by a solid vertical line. 
The fact that the red curves converge to the black far away from the evanescent region around 
corotation is decisive evidence that angular momentum is radiated away from the BL in the form 
of density waves. 

Something else to note is that the black curve in Fig. [9] switches sign close to corotation. 
This is exactly what would be expected if angular momentum were carried by waves, since the 



angular momentum density of a wave undergoes a change in sign at the corotation radius (Binney 
fc Tremaine||2008l ). 



We also point out that in Fig. [9] the flux is approximately constant between the edge of the 
evanescent region and zu ~ 1.7, but then starts to decrease more rapidly with radius. This reduction 
in flux is caused by shocking of the wave, which is illustrated in Fig. [TUJ In this Figure we show 
the profile of the density perturbation AS/S in azimuthal direction, i.e. at fixed w, before (a, at 
w — 1.5) and after (b, at vo = 1.9) breaking for simulation 3D9d at t = 75. It is obvious that in 
the latter case the wave has shocked, with clear discontinuity in the wave profile, which converges 
towards the "N-wave" shape ( Landau fc Lifshitz| [l959). Note that the analytically-expected value 



of Cs (equation [54]), automatically takes into account the dissipation of the wave due to shocking, 
through the factor of (AS) 2 . 



We now turn to the reduction of the wave flux with radius after the wave has shocked. Fig. 11 
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Fig. 10. — Azimuthal profiles of the wave before shocking (panel a) at vo = 1.5 and after shocking 
(panel b) at w — 1.9. 



shows —dCs/dw as a function of radius for simulation 3D9d at t — 75. The black curve shows the 



value of —dCs/dzu computed from the simulation data, where Cs is calculated using equation (49). 
The red curves are the analytical expectations for —dCs/dw due to wave dissipation by shocks 



calculated using equation (55). The dashed red curve uses the value of dE/dm from equation (56) 



and the solid red curve uses the value of dE/dm from equation (57). 



In measuring AH/S in order to calculate dE/dm, we find it necessary to correct for the finite 
width of the shock in the simulation. This is because the shock in our simulations extends over 
several cells and is not a true discontinuity. Thus, we extrapolate the shock density profile from 
a simulation and estimate what the value of AS/S would be for a true discontinuity. For our 
resolution, this correction increases AS/S by w 8%. Given the stiff dependence of dE/dm on 



AS/H (leading order dependence is (AH/H) 3 , see equation (56)), this correction is necessary to 
obtain good agreement between theory and simulations. 

The agreement between the measured value of dCs/dzu and the theoretically-predicted one 
indicates that angular momentum radiated away from the evanescent region into the disk is re- 
distributed by shock dissipation. This dissipation is most vigorous after the wave has completely 
shocked around w w 1.8, but the steepest parts of the wave may shock even earlier at w < 1.8, 
so the wave angular momentum starts to be redistributed to the disk fluid even before the entire 
wave has shocked. Since ftp > ft for the upper branch outside the evanescent region in the disk, 
the wave carries positive angular momentum, and dissipation of the wave spins up the disk fluid, 



pushing matter outward, see £5.4 
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Fig. 11. — Plot of —dCs/dw for simulation 3D9d at time t = 75. The dashed black line at w = 1.8 
denotes the radius at which the entire wave has shocked. 



5.1.2. Star 



Radiation from the BL is emitted not only into the disk, but into the star as well. The 
amplitude of acoustic waves traveling into the star is proportional to p -1 / 2 , so they don't shock in 
our simulations and propagate out through the inner boundary of the simulation domain. 

The inner radial "do nothing" boundary condition is not perfectly absorbing, so waves incident 
upon it are reflected with a reflection coefficient 



£R = 



<5E, 



out 



(58) 



Here, 5S ut is the amplitude of the outgoing wave (with respect to the BL), and 5Y,r is the amplitude 
of the reflected wave. As long as cr is not close to 1, the reflected waves will not fundamentally 
affect the dynamics of the instabilities in the boundary layer, since the "overreflection" mechanism 
described in Belyaev &; Rafikov (2012) only works if €r « 1. However, the angular momentum 



current due to stresses measured by equation (49) will underestimate the true current, since the 



measured current contains components due to both the true outgoing waves (with respect to the 
BL) and the artificial reflected waves from the inner boundary. Nevertheless, the flux due to only 
the outgoing component, C^out? can be estimated in terms of the measured flux, Cs, mens, as 

,meas 



C S . 



out 



1 



(59) 



This equation can be derived by considering a linear superposition of an outgoing wave and the 
reflected wave. The reflected wave has the opposite sign of Cs as the outgoing wave and an 
amplitude that is second order in e#, since angular momentum flux is a second order quantity. 



Fig. 12 1 shows the measured current (black line) calculated using equation (49), and the 



analytically-predicted current due to waves (red line) calculated using equation (53). The analytically 
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predicted current oscillates around the measured current, which is roughly constant for w < .96. 
These oscillations can be explained by the presence of a reflected wave from the inner boundary with 
amplitude cr « .3. This amplitude can be estimated from the straightforward-to-derive relation 

r l/2 _ r l/2 

5, max S'.min /nr\\ 

^ = r l/2 , r l/2 ' ( 60 ) 
U S,max ^min 

where C# 5max an d C^min are the maximum and minimum values of the oscillations in the red, 



analytically predicted flux curve in Figure 12 1. 



The reason for the oscillations in the red curve of Fig. 12 a is that equation (53) is an accurate 



estimator of the current due to waves if there are only outgoing waves present; it does not account 



for reflected waves. However, it is straightforward to modify equation (53) so that it estimates the 
current due to only the outgoing component. To do this, we assume that there is a single outgoing 
mode present with a constant value of k m that is reflected at the inner boundary into an incoming 
mode with the same value of k w and an amplitude which is a fraction cr of the amplitude of the 
out going mode. Under these assumptions, the current due to only the outgoing component, Cs out? 



can be expressed in terms of the uncorrected current given by equation (53), Cs ;Unc , as 



_ C^unc 

1 — 2cr cot 

where ^ is a free phase parameter 



Cs,out — 7^7 ; — — 2~ 5 (61) 

1 - 2cr cos(2k m zu + y;) + e z R 



The red curve in Fig. 12 3 shows C^ out as computed from equation (61), and the black curve 



shows Cs, out according to equation (59). The oscillations in the red curve have all but disappeared 
after applying the correction to compute the flux from only outgoing waves. Moreover, the red 
and black curves are seen to be roughly constant for w < .96 and to have very similar amplitudes. 
This proves that in the star, just as in the disk, angular momentum is carried away from the BL 
by waves. However, unlike in the disk, the waves do not shock as they penetrate into the star, but 
rather diminish in amplitude to conserve energy and angular momentum. 

Comparing the scales on Figs. [9] and [12) 3, we see that the angular momentum flux into the 
star is larger than the angular momentum flux into the disk by an order of magnitude. This means 
that most of the angular momentum lost by the inner part of the disk goes into spinning up the 
star and only a small fraction is transported outward into the disk. 



5.2. Lower Branch 



We turn now to angular momentum transport for the lower branch. Fig. 13 shows Cs cal- 
culated using equation (49) for simulation 3D9e at time t = 320; the same time as Fig. [6b, when 



the lower branch is clearly dominant. Panel a of Fig. [13] depicts Cs from the inner edge of the 
simulation domain up to the inner edge of the evanescent region in the disk at w — 1.58. Panel b, 
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Fig. 12. — a) Uncorrected current due to stresses inside the star computed using equations (49) 
(black line) & ( [53] ) (red line), b) Corrected current due to outgoing waves computed using equations 
d59b (black line) & d6Tl) (red line). 



which has a different scale than panel a, depicts Cs from just inside the inner edge of the evanescent 
region to the outer edge of the simulation domain. The dashed lines in panel b depict the edges of 
the evanescent region, and the dotted line is the corotation radius (ftp ^ .42). We point out that 



although the magnitude of Cs for the lower mode in Fig. 13 is smaller than the magnitude of Cs 
for the upper mode in Fig. [l2j the amplitude of either mode is not constant in time. In fact the 
magnitude of Cs inside the star may be as high for the lower mode as for the upper mode at earlier 
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times (see Fig. 14 in £5.3). 



x 10 




(a) 




(b) 

Fig. 13. — a) Cs for simulation 3D9e at time t — 320 inside the inner edge of the evanescent region 
in the disk. The dashed line marks the inner edge of the evanescent region in the disk b) Same 
as a, but Cs is shown within the evanescent region and beyond it. The evanescent region itself is 
located within the dashed lines and the dotted line is the corotation radius. 



In contrast to the upper branch, it is more difficult to estimate the current due to waves using 



equations (54) and (53) in both the star and the disk for the lower branch. In the disk, the difficulty 



comes about due to the fact that there are both outgoing and reflected shocks between the surface 



of the star and the Lindblad resonance in the disk (Fig. pp). Equation (54) only applies if there are 



exclusively outgoing waves, and it is not straightforward to modify it to account for interference 
effects from reflected waves. Without modification, the plot of the analytically predicted current, 



Cs, from equation (54) would show large spikes at shock crossings, and only in between shock 



crossings when the outgoing and reflected waves are well-separated is there reasonable agreement 
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between the measured and analytically predicted values of Cs (see Fig. 5 of Belyaev et al. ( 2012| )) 



The complication with applying equation (53) to analytically predict Cs inside the star is that 
k m (to zeroth order in M _1 ) for the lower mode (£4.2.2), which makes the current estimated 



from equation (I53j) equal to zero. Thus, in order to use equation (53) we would need to know h 



to first order in 1/M, and we have chosen not to pursue such a calculation in this work. 

One thing to notice in Fig. [13] is that Cs vanishes close to the edge of the evanescent region in 
the disk. This is to be expected if we have near perfect reflection, since Cs — at the point of re- 
flection if the amplitudes of the outgoing and reflected shocks are the same. Cs ^ throughout the 
entire region between the BL and the inner edge of the evanescent region due to shock dissipation, 
which causes the shock amplitude, and consequently \Cs\, to diminish with increasing path length 
traveled by the shock. This makes Cs < between the BL and the point of reflection, because 
the shock pattern has a negative angular momentum density inside corotation, and the outgoing 
shocks, which have a larger amplitude than the incoming ones, propagate in the +zu direction. 

We can estimate a lower bound on the reflection coefficient from the evanescent region by 



comparing the scales of Figs. 13 3, and 13 3. The absolute value of Cs (for outgoing shocks only) is 
^100 times lower outside the evanescent region than inside it. Since Cs is a second order quantity 



in 5S, the reflection coefficient as defined in equation (58) is cr > .9. This is only a lower bound on 



the reflection coefficient, since outside the evanescent region Cs could be dominated by an upper 
branch component. 

We also point out that Cs changes sign within the evanescent region close to the corotation 
radius in the disk (vjcr — 1-78), as expected for angular momentum transport by waves. Cs does 
not change sign at the corotation radius located in the BL at w « 1, since both the direction of 
wave propagation, and the sign of the angular momentum density change there. In other words, 
inside corotation at w w 1 the angular momentum density is positive, but waves propagate in the 
—w direction so Cs < 0. Just outside it, the angular momentum density is negative (since ftp < Q 
there), but waves propagate in the +£u direction so Cs < as well. 



5.3. Relation of Angular Momentum Transport to Temporal Evolution 

We now discuss how angular momentum transport in our simulations controls the temporal 



12 



to 



evolution of the system. One can use a characteristic value of \Cs\ ^ 3x 10 -3 taken from Fig. 
estimate a typical "evolution time", for the system when the upper branch is dominant. We define 
tevoi as the time it takes for the evanescent region (upper branch) in the disk to lose half its angular 
momentum. Taking the the evanescent region for upper branch to extend from 1 < w < 1.25, we 
estimate 



W>1 ~ 300. 



(62) 
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Fig. 14 shows a radius-time plot of Cs for simulation 3D9a. From the figure, it is clear that 
tevoi is within a factor of several of the span of time during which the upper branch is dominant. 



Moreover, we see as well from Fig. [14] that t ev oi is a l so a good estimator of the timespan during 
which there is a significant angular momentum current radiated away from the BL and the accretion 
rate is highest. 



2.5 



"D 

03 



1.5 



upp^ branch 



lower branch 



. i .TtT.rrxrrrrrrrrn. 

50 100 150 200 250 

time 



' i I i i i i I i i i i I i 
300 350 400 



x 10" 


-1 



Fig. 14. — Radius-time image of Cs for simulation 3D9a. The horizontal bars show the periods of 
time when the sonic instabilities are active, the upper branch is dominant, and the lower branch is 
dominant. 

The fact that the evolution time is comparable to both the timespan over which the upper 
branch is dominant and the timespan for significant mass accretion in the simulations suggests the 
following scenario for the temporal evolution of the system. The amplitude of the upper branch is 
initially high, but diminishes when it has depleted the evanescent region in the inner part of the 
disk of angular momentum. The system then switches to lower branch behavior with continuing 
mass accretion. However, without a mechanism of resupply for the inner disk, which does not exist 
in our simulations, the rate at which the inner disk is depleted of angular momentum must drop. 
Thus, the angular momentum current and the mass accretion rate both plummet for t > t evo i as 
seen in Figs. [l}) and [TIJ This occurs in step with a reduction in the amplitude of the acoustic 
modes, which are responsible for the angular momentum transport. 

In contrast to our simulations, there does exist a mechanism of resupply for the inner disk in an 
astrophysical system, namely accretion due to turbulent stresses. These stresses are not present in 
our simulation, since we do not have magnetic fields, which lead to the MRI instability in the disk. 
However, one may speculate that if the accretion time due to turbulent stresses were comparable 
to or shorter than t e vob then the acoustic modes would operate at a high amplitude indefinitely. 
Indeed, the amplitude of the modes themselves could potentially be controlled by the time for MRI 
to resupply the inner part of the disk. 
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5.4. Relation of Angular Momentum Transport to Mass Accretion Rate 



Angular momentum carried by the waves launched in the BL is eventually deposited into the 
disk fluid as a result of wave dissipation. In this section we relate the angular momentum transport 



rate to the mass accretion rate through the disk, M, defined in equation (25) 



We can use equations ([15]) and ([16]) to write M in terms of Cs as 

d 



M 



dw 



on dc s 



dt dw 



(63) 



Time evolution of ft represented by the first term on the right hand side forces gas to move radially, 
since to conserve angular momentum fluid elements must readjust their radial distance, m, when 



ft varies. The dft/dt term in equation (63) vanishes for a steady state disk, but turns out to be 
important for our simulations, since the rotation profile evolves in time (Fig. [2]). 



We plot in Fig. (15) the mass accretion rate in the BL as a function of radius for simulation 
3D9a averaged between times t = 70 — 100, when the upper mode with high ftp clearly dominates. 
The solid red and blue curves show the contributions of the dCs/dw and 8ft /dt terms, respectively, 
to the mass accretion rate (equation |63|), as measured from the simulation, and the black curve 



shows the total mass accretion rate measured directly from the simulation using equation (25). 



It is clear that both the 3ft /dt and dCs/dzu terms contribute significantly to the mass accretion 
rate, even far out in the disk. As a consistency check, the sum of the solid red and blue curves is 
shown in purple and is seen to equal the black curve up to small wiggles that are artifacts of the 
measurement procedure. 



Fig. 15 clearly shows that the dQ/dt term is important for the calculation of M. In our case, it 
is significant because the inner disk is depleted of mass as a result of angular momentum transport 
by acoustic modes. This causes time-varying density and pressure gradients in the disk, which is 
significant even quite far from the BL, at separations ~ w* from it. Time evolution of this pressure 
support ultimately causes explicit time dependence of Vl and drives accretion even in the region 
where the wave has not shocked yet (interior to w « 1.8 in Fig. 15). 



Prior to its dissipation the wave propagating through the disk cannot affect its state (Goldreich 



& Nicholson 1989). Only after it shocks, transfer of its angular momentum to the disk drives 



additional mass accretion. The dashed red curve in Fig. (15] shows the theoretically predicted value 
for the component of mass accretion due to the dCs/dzu term in equation ( |63| ), assuming the 
dCs/dw term is entirely due to shock dissipation. In other words, dCs/dw is computed using 



equation (55), where we use equation (57) to calculate dE/dm, and measure AH/S in the same 



way described in £5.1.1 We point out that the formula for the dCs/dw term due to shocks for 



the lower mode was given in Appendix B of Belyaev et al. (2012), and the derivation for the upper 
mode is nearly identical, see $5) 



The good agreement between the solid and dashed red curves for w > 1.8 in Fig. 15 shows 



that the dCs/dw component of the mass accretion rate is well-described by dissipation in shocks. 
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For smaller values of vj, the waves emitted from the BL have not yet fully shocked (j |5.1.1 ), so the 
agreement is not as good. 
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Fig. 15. — Plot of M measured from the simulations (black line) and the contributions to M from 



the dCs/dzu (solid red line) and dft/dt (blue line) terms in equation ( |63| ). The dashed red line 
shows the theoretical estimate for the dCs/dzu term assuming wave dissipation in shocks. 



Fig. 15 



It is interesting to note that the mass accretion rate due to wave dissipation is negative in 
i.e. the disk material is pushed away from the BL beyond the point where the wave has 



shocked. This is due to the fact that Fig. [15^ corresponds to a time when the upper branch is 



dominant. As discussed in §4.2. 1\ the upper branch has a fast pattern speed with a corotation 
radius in the disk at w cr « 1.12 for this simulation. Outside of corotation (zu > zu cr ), ftp > so 
the wave has a positive angular momentum density. Thus, dissipation of the wave increases the 
angular momentum of the fluid elements in the disk, causing them to move outwards and resulting 
in a negative mass accretion rate. 

The situation is different for the lower and middle branches. Unlike the upper branch which 
has an evanescent region directly adjacent to the boundary layer, the middle and lower branches 
have their evanescent regions in the disk separated from the boundary layer. Waves can propagate 
between the BL and the evanescent region in the disk. Since ftp < ft in this region, this implies 
that the wave has negative angular momentum density and dissipation of the wave results in a 
positive mass accretion rate. Beyond the evanescent region in the disk, ftp > ft for both the lower 
and middle branches so the mass accretion rate due to wave dissipation is in principle negative in 
this region, just as for the upper branch. However, very little flux penetrates the evanescent region 
(j |5.2| ), so the contribution to mass accretion from wave dissipation beyond the evanescent region 
in the disk is very small for the lower and middle branches. 
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Effect of Stratification on Acoustic Modes 



Up to now we have discussed the results of unstratified simulations, and we find that den- 
sity stratification does not produce fundamentally different behavior. In particular, the upper, 
lower, and middle branches are all observed in stratified simulations, as well, and facilitate angular 
momentum transport. 



Panels a and b of Fig. 16 show snapshots of the upper and lower branches for simulation 3D9c 
at times t — 160 and t — 320 in the color variable w^/Ya^v^) z . Comparing Fig. 16 with Fig. [6j we 
see that the spatial morphology of the wave branches is the same in the stratified and unstratified 
cases. In addition, we showed in ^4] that the dispersion relations ( [32] ) (upper branch), ( [38] ) (lower 
branch), and ( [42] ) (middle branch) provided a good estimate of the pattern speed in simulations 
irrespective of stratification. Thus, we conclude that stratification has little effect on the dynamics 
and morphology of the modes for an isothermal equation of state. 
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Fig. 16. — Panels a and b show snapshots of w\/T^{v vu ) z when the upper and lower waves are 
dominant, respectively, although there appears to be a superposition of both types of waves in the 
panels. The dashed and dotted circles depict the boundaries of evanescent regions and corotation 
radii, respectively. In panel a, the evanescent region in the disk extends all the way to the BL. 

Although stratified and unstratified simulations exhibit similar behavior, there are some ob- 
vious differences. For instance, in stratified simulations we observe the formation of a neck in the 



density profile, which is shown in Fig. 17 The formation of this neck implies that sonic instabilities 
operate preferentially near the equator and spin up an equatorial belt on the star. 



Further evidence in favor of a spun up equatorial belt can be seen in Fig. [18] which shows 
contours of the density- weighted azimuthal average of at t = (panel a) and t — 150 (panel b). 
It is clear from panel b that the equatorial region of the star is spinning faster than higher latitudes. 
For an isothermal equation of state or more generally for any equation of state that is a function of 
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Fig. 17. — Panels a and b show contours of the azimuthally- averaged density profile at t = and 
t = 150, respectively, for simulation 3D9c. Formation of a neck in the density distribution is clearly 
visible in panel b. 



the pressure alone, deviation from rotation on cylinders as see in Fig. [18p implies the existence of 
meridional circulation. Thus, although we do not see evidence of significant azimuthal spreading of 
disk material over the surface of the star, we do observe the sonic instabilities to induce meridional 
flows. 

In principle, one may also interpret the vertical extent of the equatorial bulge of the spun-up 



fluid in our stratified simulations as evidence for the emergence of a spreading layer (Inogamov & 



Sunyaev|1999 ; Piro fc Bildsten|2004a ) on the stellar surface. However, given the simplified nature of 
our runs (e.g. isothermal EOS), it difficult to unambiguously claim observation of this phenomenon 
in the simulations. 

Another difference between stratified and unstratifled simulations is that t ev oi as defined in 



equation (62) is a factor of two or so higher in the stratified case. This implies that the acoustic 



modes have lower amplitudes and take a longer time to transport angular momentum away from 
the inner disk in the stratified case. Thus, mass accretion and angular momentum transport can 
occur at an appreciable rate for a longer time before shutting off when the inner disk has been 
depleted of material (j|5.3|). 



Panels a and b of Fig. 19 show radius-time plots of (v m ) and C#, respectively, for the stratified 



simulation 3D9c. These radius-time plots can be compared with those for the unstratifled case 
shown in Figs. [TJd and 14, Both the unstratifled and stratified radius-time plots look fundamentally 
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Fig. 18. — Panels a and b show contours of the density-weighted azimuthal average of at t = 
and t = 150, respectively, for simulation 3D9c. There is a slight deviation from rotation on cylinders 
in panel b with the equatorial plane rotating faster than higher latitudes. 

similar. However the features in the stratified simulation appear to be "stretched" in time as 
compared to the unstratifled one, providing evidence that t evo l is longer in the stratified case. 



7. The Incompatibility of a Models and Angular Momentum Transport by Waves 

As we have shown, angular momentum transport in our simulations is facilitated by quasi-2D 
acoustic waves, rather than turbulent stresses. This assertion follows from the fact that Cs is 



well-described by theoretical formulas applicable for waves §5.1[ and that most of the power in 
the fluctuations is in k z = modes §4.4| Because a-models are based on the notion of a local 
effective viscosity due to turbulent stresses, whereas waves are inherently nonlocal, transporting 
angular momentum over long distances before they are absorbed or dissipate, it is clear that any 
conventional model of a viscosity is inapplicable to the BL. 

We can make this point more concrete with three specific reasons why an a model is not 
applicable to transport by waves. First, in conventional a models the stress vanishes (and Cs 
changes sign) when dfl/dvj = 0, and this is not observed in our simulations. Instead, in a wave 
model, the stress vanishes (and Cs changes sign) when = fip, and we do see evidence of this (see 



I; 5.1.1 and 5.2). Although there is a corotation resonance in the BL, which could in principle be 
located close to the point where dVt/dvo — 0, this is a fundamentally different criterion. Moreover, 
for waves there is another corotation resonance in the disk, where Cs changes sign, which is an 
effect not captured by a models. 



Second, equations (53) and (54), which describe transport by waves are fundamentally different 
in structure from equations (26)-(28) and depend on quantities such as AH/S and ftp. The former 
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Fig. 19. — Panels a and b show radius-time plots of (v w ) and Cs for simulation 3D9c. 



is a function of time that depends in a complicated way on the details of the sonic instability, and 
the latter depends on the particular wave branch. Thus, although possible, it would be contrived 
to boil down an expression containing these quantities into an effective viscosity parameter; any 
such parameter would be a complicated function of the wave branch (upper, lower, or middle), 
azimuthal wavenumber, wave amplitude, Mach number, radius, etc. 

Third, Navier Stokes viscosity with an enhanced value of the viscosity coefficient, z/, due to 
turbulence assumes local transport of angular momentum, which is the reason why the stress is 
proportional to dVt/dvo. When angular momentum is carried by waves over long distances, the 
transport is nonlocal — different parts of the disk can exchange angular momentum only when 
waves dissipate, which often happens at large separation from their launching site. Thus any model, 
such as a viscosity, which is based on Navier-Stokes viscosity no longer applies. The possibility of 



non-local momentum transport has been hinted at by Inogamov & Sunyaev (2010), who postulated 



that angular momentum transport within the spreading layer (a BL with substantial extent in the 
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6 direction) of a neutron star occurs by wave damping, in order to avoid excessive heating at the 
base of the layer. 



For these reasons, we believe it is better to use equations (53) and (54) as they are for estimating 
angular momentum transport in the BL. We do not provide an effective a for our simulations, since 
it is incorrect to adopt an a viscosity model when angular momentum is transported by waves. 



8. Discussion 

The major result of this work is a new picture of the momentum, energy, and mass transfer 
in astrophysical BLs. We have found that the isothermal boundary layer naturally supports three 
k z = branches of acoustic modes when there is a supersonic velocity drop across it. These three 
types of acoustic modes are directly related to the three types of waves found for the isothermal 



vortex sheet in Belyaev &; Rafikov (2012). Each of the three acoustic mode wave branches is 
associated with a different spatial morphology. However, two features common to all of the wave 
branches are emission of sound waves from the BL into both the star and the disk and the presence 
of two corotation radii - one within the disk and the other within the BL. The corotation radius 
within the disk lies inside a broad evanescent region (see Fig. [6]). 

The behavior of a particular mode is affected by the density stratification inside the star and 
the differential rotation inside the disk. The latter can potentially lead to the formation of a pattern 
of trapped shocks inside a resonant cavity in the disk near the BL. For this to occur, the inner edge 
of the evanescent region inside the disk must be well-separated from the BL. 

An important point is that the dispersion relations for all three acoustic mode wave branches 
remain the same in 2D, 3D unstratifled, and 3D stratified simulations. This can be seen from 
Tables [2] & [3] and Figure [5j Thus, the proposed picture of angular momentum transport in the BL 
is robust with respect to the dimensionality of the problem, strengthening the results previously 



obtained by Belyaev et al. (2012). 



An acoustic mode with pattern speed ftp couples to fluid at the corotation radius within the 
BL (fl{w c ) = ftp), sapping it of angular momentum and driving accretion onto the central object. 
Thus, the fluid in the BL loses angular momentum to waves, which propagate away from the BL. 

A fraction of the angular momentum radiated away from the BL is redeposited in the disk 
by wave dissipation in shocks, which drive the evolution of the disk surface density (j |5.4| ). An 
additional factor in the time-evolution of the disk surface density is the time-dependence of fi, 
which is connected with the temporal evolution of the radial pressure gradient in the disk. The 
time-dependence of ft is an important effect in our naturally-evolving simulations. However, in a 
truly steady-state model of the BL, this contribution to the temporal evolution of the disk surface 
density would be absent. 



The remainder of angular momentum lost by fluid in the BL is radiated into the star. Waves 
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propagating away from the BL and into the star weaken in amplitude as a result of angular momen- 
tum conservation in a medium with a rising density. For this reason, waves inside the star do not 
shock — their nonlinear evolution is too weak and in our simulations they simply propagate out 
through the inner boundary of the simulation domain. In a real system, they would presumably be 
damped by some linear dissipation process such as conduction or radiation, ultimately transferring 
their angular momentum and energy to the star. 

The radiation of angular momentum by waves away from the BL and its subsequent redis- 
tribution by damping of these waves are manifestly nonlocal processes, see fF| This casts serious 
doubt on the utility of semi-analytic models employing an (intrinsically local) (^-prescription for 
viscosity in predicting BL structure and observational manifestations. 

The general picture of acoustic wave-mediated BL dynamics advanced in this work may po- 
tentially have several important consequences that we discuss next. We warn, however, that until 
simulations with a more realistic equation of state are performed, some of the following predictions 
may be somewhat speculative. 



8.1. Implications: Energy Transport in the BL 



Waves propagating away from the BL carry not only angular momentum but also energy, which 
is released by accreting material as it is spun down. Wave damping in the disk or star results in the 
nonlocal release of this accretion energy — energy is deposited into the disk or star at the location 
where the wave damps. For example, the rate of thermal energy release in the disk per unit radial 
distance due to dissipation of waves with pattern speed ftp is 



dEd 
dvo 



[Q P -Sl(m)] 



dC s 



(64) 



where dE^/dzu denotes the energy dissipation rate per unit radius (i.e. (dEdfdzu)dw is the energy 
dissipated in an annulus between w and w + dw). By our definition, which is consistent with the 
definition in §6.4.2 of Binney & Tremaine (2008), Ed < and is non-zero only where dCs/dw ^ 0, 



(i.e. after the wave has shocked in the disk). This non-local energy transport may affect the 
structure and observational manifestations of the BL in the following way. 



Conventional wisdom based on local a-models (e.g. |Balsara et al.| ( |1999[ )) states that the 
intense energy release has to take place in the BL itself, heating the gas there and considerably 
modifying the structure of the BL. This does not need to be true in our picture of wave-mediated 
transport: kinetic energy of material accreting through the BL can be transported away from the 
layer by waves and deposited in both the disk and the star. This can easily prevent the layer from 
overheating, meaning that its structure may be affected by accretional energy dissipation weaker 
than was thought before. This is a very important consequence of wave-driven transport. 



Unfortunately, we cannot explore the effect of nonlocal energy transport on BL structure in 
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our simulations directly, since they employ an isothermal equation of state and may not accurately 
capture the details of wave dissipation, especially inside the star. However, the removal of accretion 
energy from the BL in the form of waves that is a prediction of our model adds some support to the 
results of our purely isothermal simulations, which assume the fluid temperature does not change 
at all. 



As we have seen in £5.1 , at least during some periods of time \C S \ going into star is much larger 



than that going into the disk. Since the energy density of a wave is directly proportional to its 
angular momentum density, E w = QpL w (Binney & Tremaine 2008), it is natural to expect that 



most of the accretional energy released in the BL is ultimately deposited in the star, rather than 
in the disk. As a result, heating of disk material by wave dissipation, which is distributed over a 
radial range Aw ~ w* is unlikely to exceed the energy release in the inner disk due to conventional 
anomalous viscosity, such as MRI. The latter does not operate in our simulations since we have 
focused on the pure hydro case. Thus, if acoustic waves dominate transport in the BL one might 
expect the effective temperature of the inner disk and its overall spectrum not to change much 
compared to models not accounting for the presence of a BL, which seems to be precisely what is 
observed in some systems ( |Ferland et al.|[l982 ; Polidan et al.||1990 ). 

Depending on the nature of the accreting object, deep heating occurring inside the star as 
a result of wave dissipation may affect its internal structure. This effect is likely to be minimal 
for white dwarfs and neutron stars, which are supported by degeneracy pressure, but may be 
important for accreting pre-main sequence stars. Indeed, a typical T Tauri star accreting at the 
rate M = 1O _6 M yr _1 (typical during the so-called FU Orioni outbursts) should have a rate of 
accretional energy release in the BL comparable to its own luminosity ~ Lq. Depending on the 
depth at which the waves excited in the BL dissipate, the nonlocal energy transport into the star 
may be able to affect its structure ( Prialnik fc Livio||1984 ; Hartmann et al.||1997 ; Popham||1997 ). In 



fact, Kenyon et al. (1988) and Popham et al. (1996) have previously shown that during FU Orioni 
outbursts, T Tauri stars may increase their radii by a factor ~ 2, which might be explained in our 
model by accretional energy release inside the star mediated by acoustic waves. 

Observational signatures of the BL should also be affected by nonlocal transport of energy 
in the form of waves. Accretional energy carried away by waves is either buried in the star or 
transferred to the disk, both of which affect the spectrum of the BL. Given that wave dissipation 
is expected to only mildly heat the inner disk above its standard temperature set by the visco- 



turbulent heating rate (Shakura &; Sunyaev 1973), it may even be natural that the spectrum of 



the system will have a rather weak contribution from the BL as most of the liberated energy is 
buried in the star. In the case of protostellar accretion, the sequestration of accretional energy 
released during FU Orioni outbursts inside the star may help explain the so-called "luminosity 
problem" — an apparent discrepancy between the expected accretion luminosities of protostars 
and their observed luminosities ( Kenyon et al.|[l990 ; Evans et al.||2009 ). BL formation is expected 
during FU Orioni outbursts with their high value of M, which allows the disk to overcome magnetic 
truncation. A sizeable fraction of the energy released during these intense bursts of accretion can 
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then be deposited deep inside the star by dissipation of acoustic modes. This energy would later 
leak out of the star, on a timescale which is of the order of the thermal timescale at the depth 
where the energy carried by waves is deposited. 

The same process could potentially explain why BL emission is missing from some observa- 
tions of accreting systems (Ferland et al. 1982; Polidan et al. 1990). In this framework, existing 



observations of high-energy spectral excesses from accreting systems, such as cataclysmic variables 



(Pandel et al. 1987), conventionally ascribed to BL emission, may in fact be due to some sort of 



coronal emission produced above the inner disk. Hard X-ray excesses found in accreting neutron 



star systems can be due to emission from a spreading layer (Gilfanov et al. 2003; Revnivtsev & 



Gilfanov||1984 ), which has a different morphology from the BL. 



8.2. Implications: stellar spinup 



Deposition of angular momentum carried by waves deep into the star will cause stellar spinup 
to occur in a non-trivial fashion — angular momentum will be added to the star where waves 
dissipate, spinning up the inner parts of the star, which are hidden from direct view. As a result, 
the surface layers of the star may rotate slower than would be expected if the angular momentum 
of accreted material were deposited locally, near the BL. Unfortunately, we cannot directly explore 
this issue in detail in our simulations due to their simplified nature. 



Previously Barker h Ogilvie (2010) have proposed a similar picture of inside-out stellar spinup 
by dissipation of tidally excited internal gravity waves near the stellar center. Clearly, this mode 
of stellar spinup is quite different from the predictions of standard BL models based on the notion 
of local a- viscosity (Kippenhahn fc Thomasp978). 



8.3. Implications: mixing in the BL 



The standard paradigm of local viscosity in the BL generally assumes that angular momentum 
transport in the BL is facilitated by turbulence resulting from nonlinear saturation of some local 



instability such as Kelvin-Helmholtz (Kippenhahn & Thomas 1978; Inogamov & Sunyaev 2010), 



baroclinic (Piro & Bildsten 2007), etc. Quite naturally, turbulence in the layer would also result 



in efficient mixing of fluid elements, leading to elemental homogeneity if different parts of the layer 
have different chemical compositions. 

The situation is different if transport in the layer is facilitated by acoustic modes as we find in 
this work. Wave-driven transport in the supersonic case does not need an exchange of fluid elements 
to occur and thus does not necessarily result in efficient mixing inside the layer. This agrees with 
what we generally find in our simulations, which typically exhibit large-scale regular structures like 
the pattern of trapped shocks and show very little evidence for small-scale turbulence that would 
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efficiently mix the layer. Thus, one may expect elemental mixing to be less prominent in our picture 
of the BL, compared to conventional models. 

Having said that we point out that previously |Belyaev et al. (2012) found using 2D simulations 
that the regular pattern of fluid motion in the vicinity of the layer can be occasionally interrupted 
by the short-term bursts of very irregular, rather turbulent behavior, which might result in efficient 
mixing of the layer. However, we do not observe such bursts in our 3D simulations, possibly because 
in 3D as opposed to 2D, our simulations do not develop an inflection point in the rotation profile 
within the BL. This important issue clearly needs further study. 



8.4. Implications: variability of emission 

Another implication of transport by waves is that they naturally produce temporal variability, 
since a pattern of acoustic modes rotating at a well-defined pattern speed invariably leads to periodic 
emission from the BL. Note that this mechanism of variability is completely different from previous 



suggestions such as a non-axisymmetric bulge in the disk near the BL (Popham 1999) or shallow 



surface waves on a spreading layer (Piro & Bildsten 2004b), which can provide additional sources 
of variability. 

Our simulations do not implement radiative transfer and thus cannot directly link periodicity 
of fluid patterns to spectral variability of the BL. Nevertheless, temporal behavior of the BL has 
the potential to provide the most direct link between our work and observations, since from our 
dispersion relations for the upper, middle, and lower branches it is possible to work out the fre- 
quency spectrum of the BL. It is clear that it should feature a peak at frequency uj — mfip, which 
corresponds to the apparent periodicity of the wave pattern to an external observer. 

In our simulations, we typically find quite high values of uj — mVtp ~ (4 — 6)fi(iu*) for the 



most commonly observed pattern of trapped lower branch modes. Moreover, from equation (39), 
we have the following lower bound for the frequency of the lower branch mttp > MQ(vu ic )/2. This 
lower bound makes it unlikely that acoustic modes (or at the very least the lower mode, which 
is prevalent at late times in our simulations) can explain the origin of dwarf novae oscillations 
(DNOs) in cataclysmic variables ( |Warner [2 004), which typically have M ~ 100 and mQp < 



( |Patterson |[l98T] ) . 



On the other hand, the nature of our simulations may prevent us from seeing all possible 
sources of variability in the BL. In particular, a well known peculiarity of the isothermal equation 
of state used in this work is that the Brunt- Vaisala frequency 

iV 2 = ( 7 -l)4 (65) 

is zero if 7 = 1. This means that gravity waves (g- modes) in the star are neutrally buoyant (uj = 0) 
and thus do not show up in our simulations. 
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A more realistic equation of state having 7 > 1 would be able to support both acoustic modes 
(j?- modes) and gravity waves (g- modes) inside the star. The latter correspond to the lower ("— ") 
sign in the dispersion relation ( |34| ) and have lower frequency than acoustic modes, potentially 
making them relevant for explaining DNOs. In principle, density waves in the disk could couple to 
both g- modes and p- modes in the star, whereas only the latter coupling is possible for the isothermal 
equation of state. Thus, a real BL would likely support more than just the three acoustic wave 
branches described in this work and Belyaev &; Rafikov (2012), and the dispersion relation for 
gravity modes in the star coupled to density waves in the disk remains to be worked out. 



8.5. Future prospects 



Our study employs two crucial simplifying assumptions: an isothermal EOS and lack of mag- 
netic fields. Relaxing the first assumption requires using a more realistic EOS and properly char- 
acterizing radiation transfer in the system. In addition to previously mentioned consequences — 

a better treatment of thermodynamics will have 
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e.g. appearance of new modes, see §8.1 
other implications. 

For example, damping of modes excited in the BL depends on thermodynamics: it was previ- 



ously shown by e.g. Lubow & Ogilvie (1998), Bate et al. (2002) that thermal stratification results 



in wave action channeling towards the low density regions, which may accelerate the nonlinear 
evolution of modes and cause them to damp faster. This would affect the dissipation pattern of 
waves, bringing disk regions where wave energy and momentum are deposited closer to the BL and 
making the BL more pronounced from an observational point of view. Additionally, compressional 



heating of accreted material in the BL can promote the formation of a spreading layer ( jlnogamov 
fc Sunyaev||1999 ; Piro h Bildsten 2004a) on the surface of the accreting object. 



Inclusion of magnetic field in our calculations should also have interesting consequences. First, 
it would result in the presence of MRI within the disk, which will provide a means of mass delivery 
from the outer disk to the BL region. At the moment, during the periods of lower mode dominance 
in our simulations, the fluid outside of the evanescent region in the disk remains unaffected by 
acoustic modes. As a result, mass depleted from the innermost disk does not get replenished by 
external accretion. Second, MRI would introduce turbulent motions in the disk, which may interact 
with the BL modes in a non-trivial fashion and enhance the importance of k z 7^ modes for angular 
momentum transport within the disk. Third, the number of possible wave branches the system 
could support would be even higher in the MHD case, since the disk would be able to support slow, 
fast, and Alfven waves, which could couple to p and g-modes in the star. Finally, field amplification 
by intense shear in the BL would affect the internal dynamics of the layer. We leave the detailed 
exploration of these issues to future study. 
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A. Dispersion Relation for Loosely Wound Modes Having m ^> 1 in the Disk 

We derive here the dispersion relation for loosely wound modes having m ^> 1 in an isothermal 
disk (equation [32]). We assume a two dimensional setup in the w-cj) plane and ignore stratification. 
Thus, we replace the density p by the surface density S. 

We start by explicitly writing the continuity and momentum equations in cylindrical coordi- 
nates (e.g. Binney & Tremaine (2008|)), along with the isothermal equation of state. We denote 



the radial and azimuthal velocities by u and v respectively 

as 1 d , ^ N Id 
ot w ow w ocp 



<9 S 1 9 , v. X 1 9 /v. X /A \ 

+ - — {w^u) + — — (Ev) = (Al) 



du du v du v 2 d<& 1 dP 

ot ow w ocp W OW ow 

dv ^ dv ^ v dv ^ vu 1 d<& 1 dP (A3) 

dt dw w dcj) w w dcj) YiW dcj) 

P = Zs 2 (A4) 

Next, we assume a constant equilibrium disk density, which through the equation of state 
implies a constant equilibrium disk pressure. We also ignore self-gravity so there is no perturbation 
to the potential. Specifying first order quantities with a preceding 5, the linearized momentum and 
continuity equations then become 

EdSv E s 8Su 
-^ + tt^r + -^-r + -6u + Z— = A5 
ot ocp w ocp w ow 

d5u ^d5u ^ r 1 d5P ^ , k . 
_ + fi __ 2f!Jt , + __ _ o m 

dSv n dSv MzM). n , 1 dSP 

5P-s 2 5Z = (A8) 



This preprint was prepared with the A AS IAT^X macros v5.2. 
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Taking perturbations in the form exp[i J m dzu f k vu (zu f ) + zm(< 
favor of 5S, the linearized set of equations becomes 

. _ _ x 1 . x ~ im. 
zm(SZ — SZpJ— H (1 + in)ou H 

im(il — ilp)ou — zilov H — 

. _ „ , c ft 2 , is 2 ra 5X! 
zm(SZ — Up)dv + — — H — 



fipt)] and eliminating in 






0, 



where we have defined 



and used the epicyclic frequency 



n{zu) = zukm, 



//-' = 20 ( 20 + w^-\ . 

aw ) 



(A9) 
(A10) 
(All) 

(A12) 
(A13) 



We can eliminate Sv in equations (|A9j) and (|A10() by using (|Allj) to arrive at 
i m 2 (0-0 P ) 2 



2 2 

s z m z 



zu" 



5Tj m 

E tx7 



[m 2 (Vt-Vt P ) 2 - k 2 ] Su 



K 

2ft " 

s 2 m 
zu 



- (l + m)(fi-0 P ) 
[2fti-rc(ft-ftp) 



Su 



After some algebra, the above system of equations reduces to 



m 2 (n-n P ) 



9 9 

s z m z 



zu^ 



n . n 
H o - i — o 



m z 



+ 



is 



-20(* -n)- — 



(A14) 
(A15) 

= 0. (A16) 



w 2 (Q - P ) 

We now make the approximations m 3> 1 and n/m«l, and we shall refer to the latter as the 



"loose winding" approximation. In that case, equation (A16) can immediately be simplified to 

,2 



m 2 (0-0 P ) 



2 2 

s z m z 



zu* 



+ 



IS 



zu 2 (Vt-tt P ) 



-2Q(i - n) 



UK 



0. 



(A17) 



We see that the term in the square brackets is the dispersion relation (32). Therefore, it is left 



to show that the other terms are negligible. To start, note that if we take in order of magnitude 
n rsj then the "other" terms are in order of magnitude 



1 



m(Q — ftp) 



s 2 m 2 



max 



zu* 



1 n 



m rn 



(A18) 



It is straightforward to show that regardless of whether n ^> sm/zu, k <C sm/zu, or k ~ sm/zu the 
term in equation (A18) is smaller than m 2 (£l — ftp) 2 (and hence the dominant terms in the square 



brackets in equation ( A17[ )) by a factor of at least 

1 n 



max | — , — | <C 1. 

m m 



(A19) 



- 51 - 



Thus, as long asm>l and n/m « 1, the dispersion relation is accurately given by equation 
(32). In our simulations 10 < m < 40, so the first of these conditions is typically well-satisfied. It 
is more difficult to get a handle on how well the "loose winding" condition is satisfied. However, 
the fact that Table [2] shows that there is agreement between theory and simulations at the several 
percent level is good justification for the use of equation (32). 



B. Dispersion Relation for an Isothermal Stratified Atmosphere in Cylindrical 

Geometry 

We consider the dispersion relation for an unrotating stratified atmosphere in cylindrical ge- 



ometry. Our aim is to justify our use of equation (38), which is valid for plane-parallel geometry. 



For simplicity, we assume two-dimensional perturbations in the w — (f) plane. However, we do not 
initially specify the form of the gravitational potential $(07), only specifying it when necessary to 
make the calculations tractable. 



As in Appendix [Aj we start with equations (Al)-(A4). We take small perturbations on top 
of the equilibrium state, and again denote small first order quantities with a preceding 8. The 
linearized forms of the continuity and momentum equations and the equation of state are then 
given by 



1 d5Z 1 d , s . 1 dSv 
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2j at w ovj w o(p 
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where we have defined 



h s (zu) 



l_cE 
£ dm 



(Bl) 
(B2) 

(B3) 
(B4) 

(B5) 



to be the local scale height of the atmosphere. 

We now assume perturbations in the form f(w) exp[i(m(p — cot)}. Eliminating <5£ and Sv from 



equations (B1)-(B4) we have 
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Su + i 



m 



m c 



2 \ sH_P_ 
to P 



s 2 dSP s 2 SP 
-lUJOU + — — h 




0. 



(B6) 
(B7) 



P dm h s P 

Finally, eliminating Su in terms of SP we arrive at the second order differential equation in terms 
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of SP/P 



m 



8P 



+ 1 



1 d (8P 

w dvj \ P 



+ 



d 2 



dvj 2 



SP 

~p 



0. 



(B8) 



We can gain insight by comparing the differential equation (B8) to the second order equation 



for 5P/P in the plane-parallel case for constant gravity with stratification in the x-direction. The 



latter is given by Vallis ( 2006 ) as 



^1 -k 2 ) — - — — ( — 



h,dx\P 



+ 



dx 2 



Equation (B8) reduces locally to equation (|B9j) at a given radius w in the limit 

h s (w) 



(B9) 



(BIO) 



if we make the association k y — >> m/w. In our simulations, h s /w ~ .02, so the condition (BIO) is 



well-satisfied. The local equivalence of equations (B9) and (|B8[) thus provides a justification for 



our use of an effective radius in the dispersion relation (38) 



C. Dissipation Rate in the Isothermal Limit 

We derive here a formula for the energy dissipation per unit mass due to a shock in the 
isothermal limit, 7 — >• 1. This is important in the context of our simulations, since the entropy is 
not well-defined if 7 = 1 This can be seen from the expression for the entropy per unit mass for an 
ideal gas 

= ^ -^{PP' 1 ) (CI) 

dm T 7 ( 7 -l) v H J v J 

is ill-defined if 7 = 1. 

The energy dissipation rate of the wave used in equation (55) for an isothermal shociQis given 



by 



where the A denotes the difference in a quantity across the shock. It is important to show that 
AdS/dm and hence dE/dm converges to a well-defined value when 7 — >> 1. We perform this below, 



and derive an exact result for the dissipation rate dE/dm in the isothermal limit. Larson (1990) 



1 An isothermal shock means that the temperature is the same on either side of the shock but not necessarily that 

7 = 1. 
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has given an approximate expression for dE/dm valid for arbitrary 7 in the limit of weak shocks. 



We note that our exact expression (equation [57]) matches the expression of Larson (1990) when 
7 = 1 (equation (56] ) to leading order (third) in AS/S. 



To start, we combine equations \C1\ and \C2\ to write 

dE_ s 2 
dm 



In 



(El 



(C3) 



7(7 - 1) 

where the subscripts "0" and "1" denote pre and postshock quantities respectively. Using the 
Rankine-Hugoniot relations and defining 7 = 1 + 8 and e = pi/ Pq — 1, equation (C3) becomes 

,2 r 2 + e (2 + ( 5) 



dE_ _ 

dm ~ (1 + 5)8 



In 



-(1+5) 



(C4) 



2 - e< 5 

Next, we take the isothermal limit S — > and expand the expression inside the logarithm in powers 
of S. Equation (C4) then becomes 

„2 



dE 



dm (1 + 8)5 



In 



t(2 + e) (1 + t ) h(1 + e) 

2(1 + e) ' v ' 



(C5) 



Expanding the logarithm in powers of 8 and keeping the highest order term we arrive at equation 

2 £(2 + e) -2(l + e) ln(l + e) 



(57) 



dE 
dm 



2(1 + 6) 



(C6) 



Note that this expression is valid for all values of e. 



